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Abstract 

o 
o 

■ A lattice Boltzmann model with interacting particles was developed 

• in order to simulate the magneto-rheological characteristics of mag- 

C$ ' netic fluids. In the frame of this model, 6 + 1 species of particles are 

wo: 

O ■ allowed to move across a 2D triangular lattice. Among these species, 



6 of them carry an individual magnetic dipole moment and interact 
themselves not only as a result of local collisions, as in current Lat- 



^ ■ tice Boltzmann models, but also as a result of nearest neighbours 



magnetic dipole-dipole interaction. The relative distribution of the 
individual magnetic moments is determined by the intensity of an 
external static magnetic field acting on the whole system. 

This model exhibits some relevant characteristics of real mag- 
netic fluids, i.e., anisotropic structure formation as a result of mag- 
netic field induced gas-liquid phase transition and magnetic field 
dependence of the sound velocity and the attenuation coefficient. 
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I. INTRODUCTION 



Magnetic fluids, also known as ferrofluids, are ultrastable colloidal suspensions 
of subdomain ferro - or ferrimagnetic particles - e.g., magnetite (Fe^O^) - dis- 
persed in various carrier liquids (e.g., water, petroleum, transformer oil, organic 
solvents, alcohols). These suspensions behave like quasihomogeneous strongly 
magnetizable liquids due to the presence of approx. 10 17 — 10 18 magnetic parti- 
cles in one cubic centimeter and unite the properties of magnetic materials with 
those of fluids in a rather spectacular way. 

Many experimental results suggested that colloidal particles in magnetic flu- 
ids always coagulate and form chain clusters, this process being enhanced in the 
presence of a magnetic field. The formation of chain clusters was observed with 
an electron microscope ||. The chain formation process, together with the reori- 
entation of individual particles in the presence of a magnetic field, are responsible 
for the anisotropy of the physical properties of the magnetic fluids 0. For exam- 
ple, the transversal magneto-optical effects (birefringence and linear dichroism) 
induced in thin magnetic fluid layers are well explained by the above-mentioned 
microstructural processes 0,0. The sound velocity and the acoustic attenuation 
coefficient in magnetic fluids are also depending on the angle between the sound 
propagation direction and the external magnetic field ||. 

Several theoretical works IPJ^^H an d even computer simulations flOHli 



were 



dedied to particle interactions in magnetic fluids. The Monte Carlo [10-12], as 



well as other computer simulations fl"3l , |l4lj performed up to date were limited to 
magnetic fluids considered only as magnetizable media, but the fluid behaviour 
of these materials were not considered using such techniques. 

The coming out of lattice gas models gave the possibility to introduce particle 
interactions (other than collisions) in the frame of simple models fi"5|-|T8f e.g., for 
immiscible fluids. The recent established Lattice Boltzmann methods having the 
capability to consider not only the rheological behaviour of multiple components 



Lattice gas models, whose first one was the FHP modfel [2^|, simulate the 
exact dynamical history of an integer number of particles moving on a regular 
lattice while conserving mass and momentum during their collisions. The single- 
particle equilibrium distribution function specifies the fluid density at each lattice 
site and also the velocity state, while equilibrium is determined by the particle 
collision rules. The current trend in cellular automata fluid simulations is to 
replace the lattice-gas (LG) approach by the Lattice Boltzmann (LB) methods 
[p!9| - |2Tl , |25| - p^1 . A Lattice Boltzmann automaton uses a real-number description 
for the particle distribution and is far less noisy than the LG approach. It is 
parallel in nature, due to the fact that all the information transfer is local in time 
and space, so that it is most suitable for the massively parallel computers. 

Following the general LB approach |^TJ , the following lattice Boltzmann equa- 
tions were considered for a fluid with S+l total components on a two-dimensional 
(2D) hexagonal lattice: 

<(x + e a ,t + l) -<(x,t)=^(x,t) (1) 
a = 0,1,..., S 
a = 0,1, ... ,b 

where n a a (x,t) is the single particle distribution function for the u-th component 
having the velocity directed along the vector e a (see below) and £l°(x, t) is the 
collision term. In order to simplify the computer program as much as possible, 
while retaining the relevant physical aspects, in our investigations we considered 
that all particles have the same mass, which is set equal to M a = 1 for all 
a = 0,1,..., S. 

The particles are located at the nodes of the 2D triangular lattice generated 
by the unit vectors e\ = (1,0) and e° = ( \ , ^ ) so that any particle posi- 
tion vector x is a linear combination of these generating vectors. Since all our 
investigations were performed on simple geometry 2D lattices having almost a 



For simplicity, we adopted a single relaxation time form for the collision term 
and so, 



T 



(3) 



where r is the mean collision time and n^ eq is the equilibrium distribution with a 
given functional form at site x and time t. The following form for n^ eq is adopted 
from m-H: 
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a=b 



n ° = n a and v° = — £ <e a 



(5) 



a=0 a=0 

are respectively the number density and the averaged particle velocity for the 
cr-th component at any (x,t) after collision and d < I is & constant. Moreover, 
when equilibrium is reached, one has u a (x, t) = u(x, t) for all a = 0, . . . S and 



S S b 

u(x,t) na (x,t) = J2J2^an a a (x } t) 



(6) 



cr=0 <t=0 a=l 

For sufficiently small \u a \, the above forms for the equilibrium distribution func- 
tions will be positive. 



B. Interaction potential 

In order to apply the general LB model to magnetic fluids, we consider the 
particles corresponding to a = as being the "carrier liquid" particles, while 
the "colloidal particles" are those corresponding to 1 < a < S. The colloidal 
particles carry a magnetic moment of magnitude m = 1, whose orientation is 
fixed during the simulation process. For simplicity, we adopted S = 6, so that 
the magnetic moments of each kind of colloidal particles are the vectors: 



temperature T, the probability of finding a colloidal particle with a given orienta- 
tion a = 1, . . . , S becomes proportional to the Boltzmann factor exp(— W 17 /kT), 
where W a is the potential energy of the a -orientation magnetic particles: 
W° p m°-H n 

W = m — = 'h{e a .H/H) (11) 

with h = fiQinH/ kT, po being the magnetic permittivity of the vacuum. 

From the above mentioned considerations, one can see that the parameters 
p, (p and h are the main characteristics for any computer run on the basis of our 
LB model. In order to study the structure formation, we always started our com- 
puter runs by assuming that colloidal particles are initially quasi-homogeneously 
distributed over the lattice with a small 1% random perturbation. Consequently, 
at t = 0, we always had 

■ „X , r -rr,^rX (, T Q!fld{% , j) 0.5 \ 

n°(t,j,t = 0) = p<Pr(H) f 1 + ^ J 

(a = 1,...,S) 

s 

n°(i,j,t = 0)=p- Y,n a (hJ,t = 0) (12) 

cr=l 

where < rand(i,j) < 1 are uniformly distributed random numbers, < i < 
ndx , < j < ndy and 

nit) = -a/tto (13) 

ESfexp( M om'--H/*T) 
At the beginning of each run, the magnetic fluid was always considered to be at 
rest, so that 

iT(i,j,t = 0) = (14) 

This implies 

n:(t, 3 ,t = 0) = ^^n°(t,j,t = 0),a=l,...,b 

n°(i,j,t = 0) = d n a (i,j,t = 0) (15) 

During the time evolution, magnetic colloidal particles are supposed to inter- 



which becomes simplified if one takes c = 1, as mentioned above. Since neutral 
"carrier liquid" particles do not interact at all, we always have G a5a = for any 
a = 1, ... ,6 when a = or a = 0. 

Having the interaction potential defined, the rate of net momentum change 



induced at each site is a simple generalisation of an expression in |21 



dv a (r f) &=Sa=b 

v. 1 ' = -rf{x,t) J2 E G^an c {x + e ai t)e a (19) 

at a=l a=l 

(the fact that all particles have the mass equal to 1 has been taken into account). 

The interaction process is achieved during the collision phase in the LB au- 
tomaton, i.e., during the collision time r a = r. The effect of the interaction is to 
modify the local velocities. Therefore, after the interaction is achieved, the new 
net momentum u° ew (x, t) at site x for the <r-th component becomes 

1 



(x,t) 



n a (x, t) 



x du a (xA) 
n°(x,t)u: id (x,t) + r- 



(20) 



dt 

where u° ld (x, t) is the local velocity before the interaction. In completely uniform 
equilibrium, there can be no relative flow of particles of different species since 
these are supposed to have the same mass. Consequently, the particle distribution 
functions must be locally proportional ]30| . Otherwise, different kinds of particles 
would have different temperature, which is unphysical. For this reason, the local 
velocity u^ ld (x,t) before the interaction should be choosen always as being the 
same for all a, i.e., one has 

V M {x,t) = u(x,t) (21) 

where u(x, t) was defined by Eq. (H) • 

The interaction potential does not conserve the net momentum at each site, 
as usual in current LG and LB methods. However, the total net momentum is 
conserved on the whole lattice. This can be seen from the symmetry properties 

of G a oa. 

Gaaa G^aa (22) 



4. The new equilibrium distribution functions n°(i,j,t) are computed in ac- 
cordance with Eqs.(f|) where u a (i,j,t) = u° ew (i, j,t). 

5. The Boltzmann equation Eq.(|l|) and the propagation step are now consid- 
ered in order to get the propagated distribution functions rf a {x + e a , t + 1). 

D. Conservation laws 

The lattice Boltzmann equation 

<(£ + e a ,t + l) - <(£ , t) = -± K(x, t) - rff q ) (24) 

can be rewritten after performing a Taylor's expansion [n a a = n^(x, t); summation 
from 1 to 2 over repeated greek indices is understood): 

dtK + (e a ) p d p n: + i (e a ) 7 (e fl ) fadpnl = - i « - n^ q ) (25) 
Retaining only the first-order derivatives and summing over a and a, one has 



d t £< + £d a (ea)X = (26) 

cr,a <r,a 

which, according to Eqs. (^) and (§), is just the continuity (mass) equation: 

<9 t p + V ■ (pu) = (27) 

The momentum conservation equation is obtained when Eq.(p5|) is multiplied 
by (e ) a and summed over cr and a; th e supplementary term —d a W was added 
because of the interaction potential: 

d t (pu a ) + dpU^ + P a = (4)a<' nc? - d a W (28) 

where 

pu a = J2(ea)a< (29) 
(r,a 

= V (e n )Je n )nn° (30) 



Since n^ eq is an equilibrium distribution function, one has d t n° ,eq = and so, 

nT eq = -r(e a ) 7 9X' £? (38) 

Consequently, we get: 



<j 



(39) 

Because of the second order derivatives, only n a a ,eq has a relevant contribution 
to P a and so, 

p ° = 2(FT2) v2 ^ J (40) 

The last term in the momentum conservation equation ( p8| ) is obtained after 
a series expansion: 

cW = ^ GUXK + (ea)A/](e fl ) a = 

cr,cr,a 

= E G a9a n a d^i 9 {e a ) a {Q p (41) 



Taking into account the expression (|Tq) of G a aa, one has 

( M -U = E GWea)a(e a )/3 = (42) 

a 

6c 2 36c 4 

— (m CT • m*)5 aP - D{D + 2) i(™ a ■ m & )8 aP + {m°) a {m*) p + (m^K)„] 
where the tensor (M^)^ is seen to be symmetric 

{M^) ap = (M^pc (43) 

Consequently, we have 

d a W = ^{M^^dprf = ^(M,,)^^) (44) 

cr,cr cr,cr 

III. SIMULATION OF STRUCTURE FORMATION IN MAGNETIC 



and the x component of the magnetisation 



s 



M x (i,j,t) = cos 



a=l 



2tt(<7- 1) 



£<(U,*) (46) 



a=0 



were retained. The value of the y component of the magnetisation was always 
found to be very small (M y (i, , j,t) ~ 0) during these computer runs, a fact 
which is a direct consequence of the x orientation of the magnetic field. The 
mean value of M x (i,j,t) 

i ndx ndy 

mm = < mm > = (ntfa + 1)(n<%+1) £ £ ".(U, (47) 

was found to be constant during the time evolution of the automaton, as expected. 

Figure 1 reproduces the automaton state after t — 0, 10, 100 , 500, 1000 and 
5000 iterations, respectively, for p = 0.5, = 0.20 and h = 0.5 (ndx = ndy = 
127). The white points in this figure have M x (i,j,t) > M x (t) + 0.1 • M x (t), the 
gray ones have M x (t) < M x (i,j,t) < M x (t) + 0.1 • M x (t), while the black ones 
have M x {i,j,t) < M x (t). The phase separation, i.e., the onset of thread-like 
clusters orientated along the magnetic field direction, is evident. Therefore, the 
white points in Figure 1 belong to the high magnetisation phase, while the black 
ones in the same figure belong to the low magnetisation phase. 

We have computed separately the mean values of the x component of the 
magnetisation in each phase, i.e., 



M^ h (t) = ^— ]T M x (i,j,t) (48) 

M x (i,j,t)>M x (t) 



Nhigh 



M l r(t) = -±- £ M x (i,j,t) (49) 

lyiow M x (i,j,t)<M x (t) 

where N high and N tow are the total numbers of sites belonging to the high and 
low magnetisation phase, respectively. A similar procedure was adopted also for 
the evaluation of the mean densities p hl9h (t) and p low (t) in the high density and 
the low density phases, respectively. The resulting values for the run in Figure 1 
are reproduced in Table 1. 

As the mean values for this run were p(t) = 0.5 and M x (t) = 0.02426, one 



(AlNiCo8) or type II superconductors, is well evidentiated in the frame of this 
model. 



B. Phase transitions 

Although the problem of the derivation of an equation of state for our Lattice 
Boltzmann model was not considered here, we made some attempts (computer 
experiments) in order to evidentiate those values of <j) and h at p — 0.5 which are 
characteristic for the phase diagram. The results of systematic searches performed 
up to date are displayed in Figures 3-5, where the field dependence of M x , M^ tgh 
and M l ° w after t = 5000 iterations was plotted for <p — 0.13, <p — 0.14 and 
= 0.15, respectively. The upper and lower magnetisation values are initially 
close to the mean value M x for lower values of the field parameter h. As the field 
parameter increases, the phase separation is achieved. This process is clearly seen 
as the bifurcation of the magnetisation curves in Figures 3-5. Consequently, 
when the magnetic field becomes greater than a critical value h c = h c (<fi), the 
high and low magnetisation values become clearly different. The critical field 
values corresponding to the concentrations in Figures 3 and 4 are h c (4> = 0.13) ~ 
1.1 and h c (<p = 0.14) ~ 0.6. For <ft = 0.15, the phase separate even at very low 
field intensity, 0.0 < h c (4> = 0.15) < 0.1, as one can see in Figure 5. 

The separation of magnetic phases at high values of the magnetic field in- 
tensity is a characteristic process for magnetic fluids 0J^,Q. This would be 
unconvenient for many industrial applications, e.g., high speed magnetic fluids 
rotary seals, but this process is partially overcome by the surfactant layer of the 
colloidal particles, which always introduces a supplementary repulsive interac- 
tion. Therefore, a more realistic Lattice Boltzmann model for magnetic fluids 
should take this aspect into consideration, in a similar way as in the Monte Carlo 
models already developed [[|-[l2[ . 



in the mass and momentum conservation laws and get (u' = u'), after taking 
into account that p eq is an equilibrium solution and retaining all terms up to the 
third order : 

dtp + p eg d a u' a + d a (p'u' a ) = (52) 



p eq d t u' a + d t (p%) + Ul ~ do)S aP d pP f + p e<? E CT f a {dpu' a u' 



D+2 



5 aP p eq d^u' + p^dpu'p + p eq d p dpu' a \ + 

p^dpdpu'a = (53) 

- \ (2p^d pP ' + 2 P 'dp P ')4> 2 rr ■ 

^(m°-m°)5 af3 - IJ g^ y [(m ff -m% + {iff) a (rn?) p + {rff) p {m & ) a 



The squared sound velocity (c^) 2 is, in the first order approximation, the 
coefficient of 5 a p(dpp'): 



,1 \2 



^ (1 - do) + P 



6c 2 



36c 4 



E rri 



m • m 



(54) 



D D(D + 2) 

This result is a generalisation of the squared sound velocity expression in current 
2D Lattice Boltzmann models |20|j2T}1 and incorporates also the influence of the 
magnetic field through the distribution functions f a , as well as the influence of 



the colloidal particle concentration (f> [23,31 



The sound propagation equation is obtained from the conservation equations 
( |52"D and (|5B|), substracting them after their multiplication with d t and d a , re- 
spectively, and taking into account the first-order approximation in the continuity 
equation: 

d 2 tP ' - (4) 2 vV + ^(vV) = 

p e v E rd a dp( u > a u'p) - p e ^ 2 (d a d p p') e rf(s a ,u + 



d a { P 'd pP ')^ e rr (M a ,) 



a 3 



(55) 



where 



i.e., 



p(i,j,t — 0) = p [1. + p cos(2ni / ndx) ] (58) 

Consequently, 

n°(t,j,t = 0)=p(t,j,t = 0)<f ) r(H) 

n°(i,j,t = 0)=p(i,j,t = 0)-(l-<!>) (59) 

The time evolution of the lattice automaton was registered over nu er = 5000 
iterations, for different values of the concentration and the field intensity pa- 
rameter h. The field direction was usually oriented along the x or y axes but, 
in order to study the anisotropy of sound attenuation, the general case when the 
angle between the field vector H and the x axis is 9, was also considered. 

The general behaviour of the space and time dependence of the perturbation 
(after a mediation over the y direction) 

'M = i^ -a (60) 

was found to be close to the the expression of the attenuated standing waves with 
the wavenumber k = 2n/ ndx 

p'(x,t) = p e~ at cos(kx) cos(uj s t) (61) 

i.e., 

p\i : t) = p e~ at cos(ki) cos(Lu s t) (62) 

for i — 0, . . . , ndx and t — 0, . . . , nj ter . 

In order to get the interesting quantities a and us = kcs, which are always 
accesible to experimental measurements, the x dependence was eliminated (L = 
2n/k is the lattice length): 

2 r L 

a(t) = - \ p'(x,t)dt = p e' at cos(cu s t) (63) 
L Jo 



C. Sound velocity 



The typical time evolution of the computed local density perturbation a(t) is 
reproduced in Figure 6, for = 0.20, h = 0.8 and two perpendicular orientations 
of the magnetic field (x and y). The different oscillating frequencies, due to the 
anisotropy of sound velocity, and also the different attenuation of the sound in- 
tensity are very clear. The computed corresponding Fourier spectrum reproduced 
in Figure 7 also illustrates the sound velocity anisotropy. From these figures, it 
is very clear that, for the same field intensity, the sound velocity is greater when 
propagating in the x direction. 

On the basis of these first results, a systematic exploration was made in or- 
der to see the influence of the magnetic field intensity at a fixed concentration 
(ft = 0.20. The Fourier spectra demonstrated that the sound velocity cs and 
consequently, also the angular frequency u s , increase in both x and y directions 
when increasing the value of the parameter h, i.e, when increasing the field inten- 
sity while temperature is maintained constant. As mentioned above, the velocity 
increase in the x direction is always greater than the corresponding increase in 
the y direction. 

Figure 8 shows the concentration dependence of the squared angular frequency 
u)g, which was obtained after performing computer runs with the magnetic field 
h = 0.8 in the x direction. One can see that the squared angular frequency is 
increasing when increasing concentration, a fact which also agrees qualitatively 
with the theoretical formula ([54]) ■ 

The general behaviour of the sound velocity vs. field intensity, which are re- 
tained by our Lattice Boltzmann computer experiments, i.e., the initial increase 
of the velocity, followed by saturation, agrees well with real experimental mea- 
surements e.g., those performed on water based magnetic fluids |32| . 



D. Sound attenuation 



V. CONCLUSIONS 



A lattice Boltzmann model with interacting particles was developed in or- 
der to simulate the magneto-rheological characteristics of magnetic fluids. In 
the frame of this model, 6 + 1 species of particles are allowed to move across 
a 2D triangular lattice. Among these species, 6 of them carry an individual 
magnetic dipole moment which becomes unchanged during the time evolution of 
the automaton. These particles interact themselves not only as a result of local 
collisions, as in usual Lattice Boltzmann models, but also as a result of near- 
est neighbours magnetic dipole-dipole interaction. The relative distribution of 
the individual magnetic moments is determined by the intensity of an external 
static magnetic field acting on the whole system. This model exhibits some rele- 
vant characteristics of real magnetic fluids, i.e., structure formation as a result of 
magnetic field induced gas-liquid phase transition and anisotropy of these struc- 
tures. The magnetic field induced anisotropy of sound velocity and attenuation 
in magnetic fluids is also well evidentiated in the frame of this model. 

The extension of this model in order to allow the particle system to be sub- 
jected to time variations of the applied magnetic field amplitude and/or orienta- 
tion through the introduction of a second relaxation time which may take into 
acccount the orientation relaxation of the magnetic colloidal particles, may serve 
as a basis for the analysis of hot topics related to the magneto - rheological be- 
haviour of magnetic fluids (surface instabilities, pipe flow, magnetic Benard con- 
vection, Taylor vortices formation, heat transfer), as well as an efficient approach 
to the simulation of the behaviour of some magnetic fluid industrial devices, such 
as rotary seals, dampers and inductive transducers, onset of the rotary motion of 
magnetic fluids under the action of rotating magnetic fields and phase transitions 
induced in magnetic fluids under the action of transient magnetic fields. 
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TABLE 1. 



Jiigh 



Jow 



Jlfhigh 



M l x ow 




10 
100 
500 
1000 
5000 



0.50000 
0.50000 
0.50003 
0.51433 
0.51276 
0.51298 



0.50000 
0.49999 
0.49996 
0.49145 
0.49083 
0.49083 



0.02431 
0.02427 
0.02435 
0.05967 
0.05440 
0.05394 



0.02418 
0.02423 
0.02414 
0.00147 
0.00070 
0.00124 



TABLE 2. 



Jiigh 



Jow 



M high 



M l ° w 




10 
100 
500 
1000 
5000 



0.50000 
0.50000 
0.50000 
0.500000016 
0.500000013 
0.500000017 



0.50000 
0.50000 
0.50000 
0.499999985 
0.499999986 
0.499999984 



0.0060778 
0.0060654 
0.0060635 
0.0060629 
0.0060628 
0.0060627 



0.0060475 
0.060599 
0.0060619 
0.0060624 
0.0060625 
0.0060626 



List of figure captions 



1. Dynamic evolution of the 
local magnetisation after t — 0, 10, 100, 500, 1000 and 5000 time steps, 
for p = 0.5, = 0.20 and h = 0.5. 

2. Dynamic evolution of the 
local magnetisation after t — 0, 10, 100, 500, 1000 and 5000 time steps, 
for p = 0.5, = 0.05 and h = 0.5. 

3. Field parameter (h) dependence of the mean (•), high (O) and low (□) 
magnetisation values after t = 5000 time steps, for p = 0.5 and = 0.13. 

4. Field parameter (h) dependence of the mean (•), high (O) and low (□) 
magnetisation values after t = 5000 time steps, for p = 0.5 and <fi = 0.14. 

5. Field parameter (h) dependence of the mean (•), high (O) and low (□) 
magnetisation values after t = 5000 time steps, for p = 0.5 and <p = 0.15. 

6. Typical time evolution of the computed local density perturbation a(t) for 
4> = 0.20, h = 0.8 and two orientations of the magnetic field (x and y). 

7. Computed Fourier spectrum for the two curves in Figure 6. 

8. Concentration dependence of the squared angular frequency cug, obtained 
for h = 0.8. 

9. Dependence of the attenuation coefficient a vs. the angle 9 for = 0.10 
and h = 1.0. 

10. Dependence of the attenuation coefficient a vs. the angle 9 for = 0.20 
and h = 0.5. 



